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rn ■ ABSTRACT 

o ■ 

Context. Many astrophysical problems, ranging from structure formation in cosmology to dynamics of elliptical galaxies, refer to slow pro- 

■ cesses of evolution of essentially collisionless self-gravitating systems. In order to determine the relevant quasi-equilibrium configuration at 
(-H , time t from given initial conditions, it is often argued that such slow evolution may be approximated in terms of adiabatic evolution, for the 
^ calculation of which efficient semi-analytical techniques are available. 

Aims. Here we focus on the slow process of evolution, induced by dynamical friction of a host stellar system on a minority component of 
"satellites", that we have investigated in a previous paper, to determine to what extent an adiabatic description might be applied. 
Methods. The study is realized by comparing directly N-body simulations of the stellar system evolution (in two significantly different models) 
from initial to final conditions in a controlled numerical environment. 

Results. We demonstrate that for the examined process the adiabatic description is going to provide incorrect answers, not only quantitatively, 
but also qualitatively. The two classes of models considered exhibit generally similar trends in evolution, with one exception noted in relation 
to the evolution of the total density profile. 
t ^ ■ Conclusions. This simple conclusion should be taken as a warning against the indiscriminate use of adiabatic growth prescriptions in studies 

■ of structure of galaxies. 
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1. Introduction 

For the study of essentially collisionless systems the paradigm 
of adiabatic growth has received significant attention in as- 
trophysics, especially in relation to two important contexts. 
Within cosmological scenarios of structure formation, it has 
been proposed as a natural framework to calculate the effects 
of infall of baryonic matter into the potential wells dominated 
by dark matter halos (see Blumenthal et al. 1986; Mo, Mao & 
White 1988; Kochanek & White 2001; Gnedin et al. 2004). In 
addition, it has been used as a tool to model stellar dynami- 
cal cusps in the vicinity of massive black holes at the center of 
galaxies (Young 1980; Goodman & Binney 1984; Cipollina & 
Berlin 1994; Quinlan, Hernquist & Sigurdsson 1995; Merritt & 
Quinlan 1998; van der Marel 1999). Either by semi-analytical 
techniques or by means of suitable simulations, it has thus been 
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shown that such adiabatic growth (usually justified by recalling 
that dissipative matter is able to accrete slowly towards the cen- 
ter) leads to a steepening of the underlying potential well. 

Recently, in contrast with the above picture, several in- 
vestigations have noted that processes usually interpreted in 
terms of dynamical friction, such as the capture of satellites 
and merging, tend to lead to systems with softer density profiles 
(El-Zant et al. 2001; Berlin, Liseikina & Pegoraro 2003, here- 
after Paper I; El-Zant et al. 2004; Ma & Boylan-Kolchin 2004; 
Nipoti et al. 2004), with a significant impact on the cuiTent de- 
bate about the observational counterparts to the universal halo 
density profiles found in cosmological simulations (Navarro, 
Frenk & White 1996; Moore et al. 1998; Ghigna et al. 2000; 
Navarro et al. 2004). 

In Paper I we set up a numerical laboratory for the study 
of dynamical friction by simulating the slow evolution of a 
stellar system as a result of the interaction with a system of 
satellites, with a distribution that guarantees that approximate 
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spherical symmetry is maintained. Such quasi- spherical model 
(inspired by the observation of globular cluster systems in some 
elhptical galaxies, but not necessarily meant to reproduce those 
conditions) has the advantage of being associated with a very 
smooth evolution, as desired if we wish to identify the small ef- 
fects that are characteristic of dynamical friction; in this sense, 
it extends and improves earUer models, based on the capture 
of a single satellite, used to study the mechanism of dynamical 
friction (see Bontekoe & van Albada 1987; Zaritsky & White 
1988; Bontekoe 1988; and many following papers). 

The process of dynamical friction is not expected to follow 
the adiabatic picture, because the "microscopic" mechanisms 
at its basis (scattering of orbits in a discrete picture or wave- 
particle resonance in a continuum Vlasov description) are in- 
herently non-adiabatic. In this paper we clarify this point by 
analyzing a set of controlled numerical experiments aimed at 
pointing out the differences between slow non-adiabatic and 
slow adiabatic evolution. For the general concepts and defini- 
tions relevant to the description of adiabatic processes in classi- 
cal mechanics, the reader is referred to classical books such as 
Landau & Lifchitz (1969) (see Chapter 49), Goldstein (1980) 
(see Chapter 11-7), or advanced monographs such as Northrop 
(1963). 

In Sect. 2 we describe the planned experiments and diag- 
nostics, both at macroscopic and at microscopic level. In Sect. 3 
we present the models and the code adopted in the numerical 
simulations. In Sect. 4 we describe the results, demonstrate the 
significance of non-adiabatic efl'ects in phase space, and show 
that the evolution of the underlying distribution function of a 
stellar system can be qualitatively diff'erent in non-adiabatic 
models from that found in the adiabatic case. In Sect. 5 we 
draw the main conclusions. In the Appendix we provide a the- 
oretical framework to the adopted diagnostics of "scatter plots" 
by introducing a suitable description of the evolution of the rel- 
evant distribution function in phase space. 

2. Two experiments of slow evolution within 
quasi-spherical symmetry 

No dynamical evolution is expected for a stellar system in a 
(stable) collisionless equilibrium. But in the presence of some 
colhsionality such steady state can be broken. This is the case 
for a stellar system that drags in, by dynamical friction, a heavy 
object or a system of heavy objects toward its center; while 
these "satellites" are dragged inwards, the reaction to the ef- 
fects of dynamical friction makes the stellar system evolve (see 
Paper I), although very slowly. The question that we wish to 
clarify by means of a set of dedicated experiments is to what 
extent such slow evolution difl'ers from a process in which a 
secular shift of matter toward the center of the stellar system 
occurs adiabatically. 

To study this problem, following Paper I we refer to a spher- 
ically symmetric stellar system hosting a shell of total mass 
Msheii fragmented into Nsheii heavy objects. In the continuum 
Umit, of Nsheii — > 00 at fixed shell mass, the system would be 
characterized by spherical symmetry; for finite values of Nsheii, 
the system is only approximately spherically symmetric. The 
synometry of the system and the focus on a shell are chosen so 



as to better identify the small efl'ects of slow evolution in the 
mechanisms that we wish to investigate. 

In the following subsections we describe two experiments 
designed to study the diflterences between adiabatic and non- 
adiabatic processes. 

2.1. A shell distribution of satellites dragged in by 
dynamical friction (a non-adiabatic process) 

The first experiment, identified as case A, is the same as for 
runs BS i and BS2 described in Paper I, that is a case when the 
evolution of the host galaxy proceeds as a result of the interac- 
tion with a discrete set of heavy objects falhng in because of 
dynamical friction; while such a shell of satellites falls slowly 
toward the center, the galaxy rearranges its stellar orbits and 
therefore changes its distribution function. 

The galaxy is taken to be described by a suitable equi- 
librium distribution function /. An initially smooth, spheri- 
cally symmetric shell density distribution of matter, associated 
with the potential <S>sheiiir), is then added to the galaxy. We re- 
call that, as described in Paper I (Sect. 3.2.2), the initial self- 
consistent galaxy potential (t>G(r), in the presence of the addi- 
tional potential <i>sheiii''), is difi'erent from that associated with / 
alone (i.e., in the absence of the shell), and has to be calculated 
separately to ensure that the initial state is indeed very close to 
dynamical equilibrium. The shell is then fragmented into Nsheii 
fragments (called also satellites) with initial velocities that are 
appropriate for test particles on circular orbits in the combined 
potential generated by the galaxy and the initially smooth shell 
(Ooir) + ^sheiiir)). 

In the experiment, the satellites interact directly with the 
simulation particles that represent the collisionless host galaxy, 
and thus experience dynamical friction, spirahng in toward the 
center on a time-scale T/r much longer than the dynamical time 
fj. The evolution of the whole system is expected to be slow, but 
non-adiabatic. 



2.2. A similar shell distribution of matter moved in 
adiabatically 

The second experiment, identified as case B, considers the evo- 
lution of the same collisionless host galaxy considered in case 
A, but in the presence of an external collisionless shell of mat- 
ter, which is moved (arbitrarily) toward the center of the galaxy 
mimicking the slow sinking of the fragmented shell of case A 
with the same time-scale tj,. » tci- In this case the dynamical 
friction process is not in action and the galaxy evolves because 
the external field experienced is (slowly) time-dependent. 

In contrast with case A, here the shell is always kept 
smooth, spherically syimnetric, not fragmented, and made to 
act on the galaxy uniquely as a slowly varying "external field". 
The evolution of the galaxy is thus expected to be, by construc- 
tion, slow and adiabatic. 
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2.3. Control cases 

We will also perform two additional experiments devised to test 
the dynamical evolution of the system when the mechanisms 
described above are not applicable. 

The first control experiment, identified as case C, is a sim- 
ulation that follows the prescription of case B, but with a much 
shorter time-scale t/o,,, so as to break, on purpose, the con- 
dition of adiabatic evolution. This experiment is meant to bet- 
ter identify those properties that are characteristic of adiabatic 
evolution with respect to those associated with the driving by a 
general external forcing. 

The second control experiment, case D, is a simulation to 
test the quality of the initial equilibrium in each of the above 
three cases A, B, and C. We consider two possibilities. One, 
in which the shell is treated as smooth, spherically symmetric, 
and coUisionless, is better suited as a control test for case A; the 
shell is here represented by a large number of simulation par- 
ticles placed on circular orbits but otherwise indistinguishable 
from the simulation particles used to represent the host galaxy. 
The other, in which the shell only acts as an "external field", is 
the more natural control test for cases B and C. Case D, in addi- 
tion to testing that the adopted initial conditions for the model 
of the galaxy and of the shell are correct, checks the accuracy 
of our simulations over long time-scales. 

2.4. Diagnostics 

To test the similarities and the differences in these cases, we 
resort to the following diagnostics. To study the problem at the 
macroscopic level, we show for the various cases the evolu- 
tion of the galaxy density profile and of the global pressure 
anisotropy, defined as 2KJKj, where and Kj- are the total 
kinetic energy of the galaxy in the radial and tangential direc- 
tion respectively. At the microscopic level, to study the detailed 
conservations in phase-space, we produce "scatter plots" (final 
vs. initial quantities) for single-particle energy E, radial action 
/, and angular momentum J. 

Ideally, in simulations belonging to case A, the three quan- 
tities, E, I, and J, are not conserved; in case B, the single- 
particle energy E is not conserved, but / and J are; in case C, 
only J is conserved. In case D, all quantities, E, I, and J, are 
conserved, and the relevant scatter plots thus define the limits 
of the numerical simulations (see Tab.Q. 

The use of scatter plots raises interesting and important the- 
oretical issues. Consider the scatter plot for the single-particle 
energy, which we already presented in the bottom right frame 
of Fig. 8 in Paper I and we will give in the upper frames of 
Fig. |3]and Fig. |4]of the present paper. This plot correlates the 
value of the energy of a single simulation particle at an ini- 
tial time til, with that of the same simulation particle at a final 
time tfi„. In an ideal code the simulation particles would fol- 
low the characteristics dictated by the coUisionless Boltzmann 
(Vlasov) equation. Under conditions for which the energy is 
conserved, any initial distribution function would originate a 
distribution of simulation particles for which the scatter plot is 
a segment, for any value of f/,„; for a finite number of parti- 
cles, the points would line up on such segment, covering the 
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Table 1. Expected conservation of single-particle quantities for 
the four cases considered (C: quantity conserved, NC: quantity 
not conserved). 



Case 


A 


B 


c 


D 


E 


NC 


NC 


NC 


C 


J 


NC 


C 


c 


C 


I 


NC 


c 


NC 


c 



range from the minimum to the maximum value of particle en- 
ergy, consistent with the adopted distribution function. Under 
conditions for which energy is not conserved, we may expect 
a distortion of that segment and, in addition, a broadening of 
the segment into a distribution of points, corresponding to a 
spread of eff'ects, at given initial energy, associated with the dif- 
ferent values of the angular momentum. Some conditions that 
determine the energy non-conservation may even determine a 
change in the value of the energy that depends on phase vari- 
ables that define the orbit of the simulation particles, at given 
initial values of energy and angular momentum. 

In real simulations, all the eff'ects mentioned above are 
mixed with numerical effects, which we may generically call 
numerical noise. In other words, undesired numerical noise 
may be responsible for distortions of the segment-type plot and, 
more naturally, a general broadening with respect to the ideal 
case. 

These considerations prompted us to move beyond the in- 
tuitive stage at which scatter plots are naturally appreciated, 
to see whether the information contained there can be set in 
a more satisfactory quantitative theoretical framework. An out- 
line of such theoretical framework is provided in the Appendix. 

3. The models and the code 

3.1. The basic equilibrium of the host galaxy 

The basic model for the initial conditions of the host galaxy is 
constructed from two different distribution functions. The first 
is an n = 3 polytrope with 

f(r,v)^A(Eo-Er-''^ (1) 

for E < Eo and zero otherwise. The single particle total energy 
is £■ = v^/2 + <i>(r) andEo = (^(R) = -G{Mc + Mmi)/R, where 
Mc and R are respectively the mass and the radius of the galaxy 
and M shell is the mass of the shell. This is the galaxy model used 
in Paper 1; it is convenient and is best suited for comparison 
with past investigations of the process of dynamical friction 
(for example, it was used by Bontekoe & van Albada (1987) in 
their pioneering work on the subject). 

The second is the f^^^ distribution function 

y(v) ^ j^^-laE+d(J^/\Ef'^-y'^] ^j) 

for £■ < and zero otherwise. Here E and J are the single-star 
energy and angular momentum. The four constants A, a, d, and 
V are real and positive: from these one can identify two dimen- 
sional scales and two dimensionless parameters, such as v and 
the dimensionless central potential ^' = -aO(Q). In this paper 
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we will focus on the case (v, ^) = (3/4, 5) (see Trenti & Bertin 
(2005) for details about this distribution function and a discus- 
sion of the properties of the resulting self-consistent models for 
different values of v and of ^'). The resulting models are more 
concentrated than the « = 3 polytrope and are characterized by 
an anisotropic pressure tensor. In general, these models provide 
a more realistic description of bright elliptical galaxies. 

3.2. The modified self-consistent equilibrium in the 
presence of a shell 

The initial shell is treated as smooth, spherically symmetric, 
and static, characterized by the density distribution; 

Psheii = A) exp [-4(r - r,HM(0)f /Rmi^ (3) 

for \r - rsheii{0)\ < Rsbeii (and vanishing otherwise); rsheiiiQ) is 
the initial radial position of the shell, where its density is pQ, 
and Rsheii is its half-thickness. 

As described in Subsection 3.2.2. of Paper 1, in the presence 
of an external potential the distribution function of the galaxy 
differs, from that in absence of it, in the potential term that now 
is the total potential of the system: <[)(r) = ^cif) + ^siwuir). 
As a consequence, it is necessary to re-compute the equilib- 
rium models of the galaxy by solving the Poisson equation for 
the self-consistent ^cir) in the presence of the external shell 
density distribution. 

3.3. The code 

The galaxy evolution is described by a coUisionless (mean- 
field) particle-mesh code, which is an improved version of the 
original van Albada (1982) code. The new code (see Trenti et 
al. (2005) and Trenti (2005)) is based on an improved Poisson 
solver scheme, in which the angular grid is dropped, while the 
radial grid is adaptive so that each cell is populated by approx- 
imately the same number of particles. 

The galaxy is sampled with simulation particles, with 
initial positions and velocities extracted from the distribution 
functions described in subsection l3.2l bv means of Monte Carlo 
methods (see Paper I). The same method is applied to derive 
the initial positions of the Nj fragments of the shell, of mass 
Mf - M shell IN f. These fragments or satellites are described by 
an extended Plummer density profile: 

GMf 



where 7?/ is the radius of each shell fragment. In the simulations 
of case A (see subsection 12. 1> . the equations that describe the 
interactions among galaxy simulation particles and shell frag- 
ments are those described in Paper 1. Each particle of the galaxy 
feels the action of the whole galaxy via its mean field and that 
of each shell fragment directly; each shell fragment feels the 
sum of the direct actions from all the galaxy simulation parti- 
cles and of the direct actions from the other shell fragments. 

In the simulations of cases B and C, there are no discrete 
objects (satellites) involved. Therefore, the simulation particles 



follow the orbits dictated by the mean field generated by the 
galaxy, as in case A, together with the force of the shell which 
is treated as an external, spherically symmetric, and smooth 
force. 

For the diagnostics mentioned in subsection l2.4l the radial 
action / of a single simulation particle, at time f, is computed 
in a semi-analytical way, starting from the values of E and J. 
The total potential <b{r) - <l>G(r) + ^sheii(r), felt by the particle, 
is first measured in the simulation at time t. Then, the effec- 
tive potential ^eff(nJ) = <!'('") + -/^/(2r^), the radial veloc- 
ity Vrin E, J) = ^J2[E - (^cff(r)], and the two turning points 
r,„in{E, J) and r,„ax{E, J) are derived. Finally, the radial action is 
found by numerical integration: I(E, J) - 2 "'" Vr{r, E, J)dr. 

3.4. Units 

The units adopted are 10 kpc for length, IO^'M© for mass, and 
10^ yr for time. Thus, velocities are measured in units of 97.8 
km/s and the value of the gravitational constant G is ~ 4.497. 

4. Numerical simulations 

The galaxy is sampled with - 250000 simulation particles 
and has a mass Mc - 2 and a half-mass radius - 0.3 (in 
code units); simulations with higher values of have also been 
run. We should recall that the issues of the reliability (stability) 
of the code and of the adequacy of the number of particles used 
in the simulations were already addressed in previous papers. 
In particular, in Paper I (see Sect. 5.1 of that paper) we ran a 
simulation of a single sinking satellite with up to 2 x 10^ and 
found very little change in the observed evolution, down to the 
innermost radii reached by the infalling satellite; in addition 
(see Sect. 4.4.2 of that paper), we checked, with general quali- 
tative agreement, the main results obtained earlier by Bontekoe 
& van Albada (1987), who had originally investigated the pro- 
cess of dynamical friction with only a few thousand simula- 
tion particles. Furthermore, the new version of the code used in 
this paper has been extensively tested with many simulations 
with A^ > 8 X 10^, especially designed to check its stability 
in reproducing concentrated models over long time scales; the 
new version of the code has also been tested against the per- 
formance of high-quality tree codes currently available (Trenti 
2005; Trenti et al. 2005). To be sure, if the process of dynam- 
ical friction depended too sensitively on detailed resonant ef- 
fects, even with these large numbers of simulation particles the 
phase space might be insufficiently well sampled. 

The mass of the shell is one tenth of that of the galaxy 
(Msheii = 0.2). The shell is initially located at rsheiiiO) - r^- 
The radius of each fragment is equal to the half thickness of 
the shell, R/ - R shell = 0.1, that is one third of the half-mass 
radius of the galaxy. These values of Msheii/Mc and of Rj/tm 
repeat the conditions of previous investigations (Bontekoe & 
van Albada 1987; Paper I); this choice helps to better bring out 
the effects of dynamical friction within the computational re- 
sources used. 

With the above choice of and Mc, the dynamical time of 
the galaxy without the shell, defined as = GM^~ /(-2E,oi)^^^, 
is f^''°'* = 0.1476 for the polytropic model and 4"' =0.1881 for 
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Table 2. The main set of numerical simulations; « = 3 for the 
polytropic model and (v, *P) - (3/4,5) for the /^''' model (see 
text for description); in runs Dl and D3, Nf represents the num- 
ber of simulation particles treated as collisionless. 



Case 


Process 


Run 


Shell 
type 


Galaxy 
model 


A 


Dynamical 


Al 


Nf = 20 


Polytrope 




friction 


A2 


Nf = 100 


Polytrope 






A3 


Nf = 20 


y(v) 






A4 


Nf = 100 


fM 


B 


External 


Bl 


External slow 


Polytrope 




adiabatic 


B2 


External slow 




C 


External 


CI 


External fast 


Polytrope 




non 


C2 


External fast 






adiabatic 








D 


Control 


Dl 


Nf = 25000 


Polytrope 






D2 


External static 


Polytrope 






D3 


Nf = 25000 


/'"' 






D4 


External static 





the f^^^ model (in code units). To better illustrate the meaning 
of td, we note that the period of two selected circular orbits, at 
radii 1.17 and 3.3 r^, is 5.4 and 14.3 for the polytropic 
model and 3.0 and 12.7 for the /'^''* model. 

In all simulations the conservation of the total energy and 
total angular momentum is respectively of the order of 10"^ 
and 10"^ per dynamical time, except for case D, for the poly- 
tropic model, in which the conservation is better by one order 
of magnitude. 

The evolution of the system is followed up to a final time, 
tfin - 9Qtd (for the fast case C, the evolution is followed up to 
tfi„ - 2td). The configuration of the system at such final time is 
compared with that at f,„ = 0.14 fj. 

The list of experiments is given in Tab.|2 

4.1. Case A: Fall of a shell of satellites by dynamical 
friction 

We have performed four simulations of this type: two with the 
host galaxy modeled by a polytropic distribution function and 
the shell sampled with 20 and 100 fragments (runs Al and 
A2, respectively) and two with the galaxy described by the /*^^ 
model and the same configurations for the shell (runs A3 and 
A4). Runs Al and A2 basically repeat runs BS i and BS2 de- 
scribed in Paper I, but with the new improved code (Paper I 
made use of the original van Albada (1982) code). 

In these runs, because of the dynamical friction experienced 
by the individual fragments, the shell of fragments sinks in 
slowly, toward the center of the galaxy, while becoming thicker 
in radius, exchanging part of its energy with the galaxy (see top 
panel of Fig.^for run Al). For the same galaxy model, by re- 
ducing the number of fragments (from 100 to 20) the fall of the 
shell is faster, as expected because the fragments have larger 
masses, and the thickening of the shell is less pronounced. 



The host galaxy responds with an expansion that leads to 
a density profile broader in the outer parts and shallower in 
the inner regions. The degree of variation of the density profile 
from the initial to the final configuration increases when the 
number of shell fragments is decreased (from about 1 1 % of 
run A 1 to about 2% of run A2, for the polytropic model) and the 
size of the effects changes with the galaxy model (from about 
ll%ofrunAl to about 25% of run A3 for A?/ = 20). The initial 
and final density profiles observed in run A 1 are illustrated in 
the central panel of Fig. [fl the plot makes use of the volume- 
weighted density r^p(r) so as to better bring out the shift of 
mass associated with the observed evolution. 

In addition, the host galaxy responds by changing its pres- 
sure tensor in the tangential direction; the effect is more marked 
in cases Al and A3, which have a smaller number of frag- 
ments. For the polytropic model, the final values of the global 
anisotropy parameter IKr/Kj (twice the ratio of the total ki- 
netic energy in the radial direction to that in the tangential di- 
rections) are 0.96 and 0.98 respectively for runs Al and A2 
(see bottom panel of Fig.^for run Al), while the initial con- 
figuration starts as isotropic (with 2K, IKt - 1). In turn, the 
model, which is characterized by a radially biased pres- 
sure anisotropy to begin with, slowly evolves toward a more 
isotropic configuration. 

Because collisions are at the basis of the underlying physi- 
cal mechanism of dynamical friction, we expect to see their ef- 
fects on the particles of the galaxy that have taken part in the ex- 
change of energy and angular momentum with the fragments of 
the shell. For the polytropic model, by comparing the first col- 
umn of Fig. |3] with the first column of Fig. 0] which represents 
the control case, one indeed sees that the single-particle energy 
E, angular momentum J, and radial action / (of the galaxy sim- 
ulation particles) are not conserved. The percentage of particles 
with quantities conserved to a given level increases with Nf. 
In the polytropic galaxy, about 66% and 45% (respectively for 
Nf = 100 and 20) of the particles have !(£■/,„ -E/J/E,-,, I < 0.08, 
where £,„ and Efi„ are the initial and final energy values, in 
contrast with the value of 94% observed in the collisionless 
control case D. The corresponding values for J and / and for 
the runs related to the ' model are shown in Tab.|3l [The gen- 
erally worse conservation of J and / with respect to that of E 
depends on our definition of conservation in terms of relative 
variations and on the fact that a significant number of particles 
have values of J and / near zero. For example, in the control 
run Dl, if we focus only on particles with J/Jmax and I/Imax 
greater than 0.25 (where J^ajc and /,„av are the maximum val- 
ues of J and / taken by the galaxy particles), the percentage of 
conservation increases respectively from 61 to 76 and from 38 
to 65 (the fraction of particles left out is 32% and 70%).] 

4.2. Case B: Fall of a shell described by a slow 
external potential 

We have performed two simulations of this type: run Bl for the 
polytropic galaxy model and run B2 for the /^''' model. 

The construction of the external force mimicking the fall 
of a shell is based on the results for the shell density distri- 
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bution with Nj - 20 observed in runs Al and A3. The exter- 
nal specific radial force to which a galaxy particle, at radius r 
and at time t, is subject is given by /r(r, f) = -GMsheiiir, t)/r^, 
where M^iir, t) = M,/,,//(r, t = 0)(f/,„ - ?)/?/,„ + M,heu{r, t = 
t fin){t / 1 fin) is the mass of the shell inside the sphere of radius r 
at time f. The values of Msheii{r, t -Q) and Msheii(f, t = f//„) are 
computed in the runs Al and A3 (for run Al, see top panel of 

Fig.ia. 

The observed evolution of the galaxy in the presence of the 
slow change of the potential associated with the external shell 
is very different from that in response to the direct interaction 
with the fragments of the shell. This is a purely coUisionless 
case in which the galaxy energy changes slowly, much like in 
cases A, but in the absence of dynamical friction. 

The density distribution becomes more peaked in the in- 
ner regions and shallower in the outer parts (see central panel 
of Fig.|2li. For a given galaxy model, the effect is found to be 
stronger when a comparison is made with a case of dynamical 
friction characterized by the lower value of Nf (Nf - 20). At 
fixed value of Nf, the effect is stronger for the /^''' model. 

The total kinetic energies of the galaxy and Kj (in the 
radial and in the tangential directions) increase with time, be- 
cause of the energy change associated with the time-dependent 
process, but at different rates, so that the galaxy global pressure 
anisotropy evolves in the tangential direction. Much like in runs 
of case A, the final configuration of the polytrope is tangentially 
anisotropic (with a final global parameter 2K, /Kt = 0.96; see 
bottom panel of Fig. |2ji. A similar trend is found for the 
model. 

At the microscopic level, as expected, the single-particle 
energy is not conserved because the potential in which the par- 
ticles move is time-dependent. On the other hand, the angu- 
lar momentum and the radial action are conserved, within the 
numerical limits (see second column of Fig. |3 and Fig. |3 and 
Tab.|3j. A similar behaviour is found for the model. 

4.3. Case C: Fall of a shell described by a fast 
external potential 

As for case B, we have performed two simulations of this type: 
run CI for the polytropic galaxy model and run C2 for the 
model. 

In these simulations the external shell is moved in, from the 
initial to the final state, as in case B, but at a much faster speed, 
that is in only two dynamical times. This rapid variation leads 
to a final configuration characterized by a density distribution 
(see central panel of Fig.|2} and global pressure anisotropy (see 
bottom panel of Fig.|5} similar to those observed in case B. In 
particular, the curves labelled Bl and CI in the middle frame of 
Fig.|2]overlap, showing that the density response of the galaxy 
to the externally imposed field mimicking the infall of a shell 
does not depend on whether the process is actually slow or fast. 
A similar behaviour is observed in the study of the model. 

At the microscopic level, for both galaxy models, the 
single-particle energy and radial action are not conserved, 
while the angular momentum is because of the spherical sym- 
metry. The numerical scatter observed in the relevant scatter 



Table 3. Conservation of single-particle quantities. For the var- 
ious runs of cases A and B the Table lists the percentage of 
galaxy particles with |(2/,„ - Qm)IQin\ < 0.08 for the quan- 
tities Q - E, J, and /. For case C the "conservation level" is 
lowered to 5%, because the scatter is smaller due to the short 
duration of the simulations. For each case (A, B, and C) the 
numbers are listed close to those of the corresponding control 
case D (cf. the expectations of Tab. The worse conservation 
of J and / with respect to that of E depends on the fact that 
a significant number of particles have initial values of J and / 
near zero. 
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D3 
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20 
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89 


I 


46 


18 


31 


48 


48 


83 


54 



plots is relatively small because of the short duration of the 
simulation (see third column of Fig.|3and Fig.l^and Tab.O. 

4.4. Case D: Control runs with a coUisionless and with 
a static external shell 

We have performed four simulations of this type: runs Dl and 
D2 for the polytropic galaxy model and runs D3 and D4 for the 
model. 

These simulations have been performed to test the stabil- 
ity of the initial conditions of the system considered in cases 
A, B, and C. In runs Dl and D3, designed to test runs A1-A4, 
the shell is sampled with 25000 particles, so that each shell 
simulation particle has the same mass as each galaxy simula- 
tion particle; the shell particles are then treated as coUisionless, 
interacting with one another and with the the particles of the 
galaxy via the mean field. In runs D2 and D4, designed to test 
runs Bl, B2, CI, and C2, the shell is treated as an external field, 
much like in cases B and C, but kept static, in its initial config- 
uration Msheiiir, t - 0). 

The initial "equilibrium" configuration of the host galaxy 
in the presence of the coUisionless or static shell is checked to 
stay basically unchanged in the course of the simulations. For 
both galaxy models, in run Dl and D3, the radii of the spheres 
containing 15, 25, 35, 45, 55, 65, 75, 85, 95 % of the total mass 
of the shell remain constant in time to within 1%, except for 
that containing 5% that remains constant only to 4%. A similar 
diagnostics for the galaxy simulation particles shows a conser- 
vation better than 0.5%; in the innermost regions (r < O.lr^), 
because of the relatively small number of particles that are in- 
volved, the density profiles are found to be steady to within 
2%. The changes in the total kinetic energy of the galaxy in the 
radial and tangential directions are found to be less than 0.3%. 
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Fig. 3. Coirelation between single-particle quantities at an initial time (f,-,, = 0.14 td, x-axis) and at a final time (f/,„, y-axis): 
energy (top row), angular momentum (central row), and radial action (bottom row), expressed in code units (see Sect. 13. 4> . The 
first column displays the results for case A (run A 1 ; f/,„ = 90 f^/), the second for case B (run Bl ; ?/,„ = 90 f^), the third for case C 
(run CI; f/,„ = 2 f^); the corresponding control cases are displayed in Fig.|3 For each quantity Q, the grey scale used represents 
the number density of simulation particles in the plane Q^-Q. By comparing each panel with the corresponding panel in Fig.|4j 
the conservation properties listed in Tab.n(and quantified in Tab.O, are confirmed. 



In runs D2 and D4 the overall stability is checked to be even 
better than for runs D\ and Z)3: the mass radii remain constant 
to within 0.05% (except for the 5% mass radius for which the 
stability is to 0.2%), the density profile remains unchanged to 
within 1%, and the total kinetic energies are stable to less than 
0.1%. 

Because of these control cases are designed to be spher- 
ically symmetric and time-independent, the single particle 
quantities (energy £, angular momentum 7, and radial action 7) 
should be conserved. The level of conservation observed over 
the entire period of the simulation, thus illustrating the accu- 
racy of our code, is shown in Fig.l^and in Tab.|3] 

For completeness, we have also run one additional control 
case, ZJ5, of a pure polytropic model, with no shell, to test 
whether some of the scatter present in the control cases D might 
be associated with the introduction of the shell itself. The rel- 
evant conservation levels over the time interval 0.14 td - 90 td. 



corresponding to the D2 column at the center of Tab.|3j are 97 
(for E), 66 (for 7), and 41 (for 7), practically identical to the 
corresponding values for run D2. 

4.5. One qualitative difference in beliavior between tiie 
polytropic and the /"^ model 

We have noted that, as a result of the slow evolution associated 
with the process of dynamical friction (see middle frame of 
Fig.Q, the density distribution of the galaxy is softened, with 
the mass being basically transported outwards, in contrast with 
the trend associated with the corresponding adiabatic case (see 
middle frame of Fig.|5}. This general property is found in both 
the polytropic and the /^''' models that we have investigated. 

If we focus instead on the evolution of the total density pro- 
file (galaxy -i- shell), we note a curious qualitative difference in 
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Fig. 4. Correlations observed in the control runs Dl and D2, defining the accuracy of the code. The layout and the notation are 
similar to those of Fig.|3] The apparent thinness of the scatter plots in the third column of this figure and in Fig.|3]are due to the 
shortness of the run presented (f/,„ - 2 td). 



the behavior of the two models. In fact, while for the poly tropic 
model (characterized initially by a broad core) the evolution in- 
duced by dynamical friction leads to a higher concentration in 
the total density profile, for the model (which starts with 
higher concentration) the softening and broadening of the den- 
sity profile occurs not only for the galaxy density distribution 
but also for the total density profile (see Fig.|5}- In all the cases 
studied, the adiabatic evolution is associated with the most con- 
centrated final density profiles. 

These considerations are especially relevant when one 
imagines astrophysical applications based on the study of rota- 
tion curves in the innermost regions of galaxies. We do not have 
a simple physical argument to explain the observed difference 
in behavior, which was noted as a side-result in this project 
dedicated to studying the non-adiabatic character of dynamical 
friction. We are planning to return to this issue in a separate pa- 
per that will be devoted to comparing the process of dynamical 
friction in models characterized by a broad core to the process 
occurring in concentrated models. 



5. Conclusions and perspectives 

In this paper we have studied in detail the slow evolution of 
quasi-coUisionless self-gravitating systems by comparing the 
situation when such evolution is due to effects of dynamical 
friction to a situation that is strictly adiabatic but otherwise sim- 
ilar The two types of evolution have been studied by specially 
designed numerical experiments and have been compared both 
at the macroscopic and at the microscopic level. 

At the microscopic level we have checked that the relevant 
single-particle quantities, such as energy E, angular momen- 
tum J, and radial action /, are indeed conserved or not con- 
served as expected from first principles. These studies have 
been quantified by means of suitable "scatter plots" which 
bring out interesting structural properties of the collision- 
less Boltzmann (Vlasov) equation, when suitably expressed in 
terms of evolution variables (see Appendix). 

At the macroscopic level, we find a qualitatively different 
behaviour in the evolution of the density distribution of the 
stellar system under investigation. In particular, in the (non- 
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adiabatic) case of slow evolution resulting from dynamical fric- 
tion, the density profile becomes "softer" in the course of evo- 
lution, in contrast with the overall steepening and contraction 
noted in the adiabatic case. In turn, the overall evolution of 
pressure anisotropy turns out to be similar in the two cases. 

In this study we have considered the evolution induced on 
the host galaxy by a slow infall of matter distributed in a rel- 
atively thin shell. This paper is part of a longer-term project, 
currently in progress, where we aim at studying the more gen- 
eral problem in which the infalling matter has a full three- 
dimensional distribution. By comparing such more general be- 
havior with the behavior that may be anticipated by superpos- 
ing the effects of adjacent thin shells of matter (now described 
quantitatively in this paper), we will be able to assess to what 
extent the general problem of infall can be reduced to a local 
analysis or else whether full global simulations are generally 
required, with a behavior that may depend on a case by case 
basis. It would be important to identify general trends in such 
general process, even if those trends will be different from those 
suggested by the picture of adiabatic infall. 

The astrophysical applications of this long-term project are 
basically directed on two goals. Specifically in the galactic con- 
text, we wish to develop a quantitative picture of the effects on 
a stellar system (such as a bulge or an elliptical galaxy) associ- 
ated with the slow infall of a system of globular clusters; sim- 
ilar effects may operate in determining a more realistic picture 
of galaxy formation, in which the capture of a large number of 
satellites (minor mergers) may accompany or follow the pro- 
cess of incomplete violent relaxation (known to occur in purely 
coUisionless systems). In the more general cosmological con- 
text, we wish to assess whether dynamical friction can change 
significantly the scenario of cusp formation, so prominently 
present in cosmological simulations. For addressing properly 
all these issues we will have to extend the calculations so as to 
allow for the evaporation and possible tidal disruption of the 
satellites that are dragged in by dynamical friction. 
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viding us with his improved code for the simulations and Luca Ciotti, 
Michele Trenti, and Tjeerd van Albada for a number of useful sugges- 
tions. We are also glad to thank Alessandro Romeo for his constructive 
criticism. Part of this work was supported by the Italian MIUR. 

Appendix A: 

We start by examining the origin of distortions and dispersion 
in scatter plots for the ideal case in which the simulation par- 
ticles follow rigorously and exactly the characteristics of the 
coUisionless Boltzmann equation. To set a quantitative frame- 
work for the corresponding effects associated with the numer- 
ical noise is far more difficult and beyond the scope of the 
present paper. The main idea that we follow below consists in 
reading the scatter plots as suitable cuts (or projections) of a 
distribution function expressed in terms of evolution variables. 

A formal definition of the scatter plots mentioned above 
can be given in terms of the distribution function / - f(x, p, f) 
of a stellar system (i.e., more precisely of its mass density dis- 
tribution function). Here x, p are the 6 position and momen- 
tum coordinates in phase space and /(x, p, f) dx dp is the mass 



inside the phase volume dx dp. The distribution / obeys the 
coUisionless Boltzmann equation 



(A.l) 



where H = H(x,p,t) is the single-star Hamiltonian that in- 
cludes the effect of external forces such as those arising from 



the presence of the falling satellites and 



lx,p 



are the standard 



Poisson brackets (referred to the variables x and p). As is well 
known, Eq. iA.l\ can be rewritten in the Lagrangian form 



Df 

— = 0, 
Dt 



(A.2) 



which expresses the fact that the distribution function remains 
constant as it evolves in time along the characteristics of the 
Hamiltonian H. In turn, the general solution of Eq. (IA.2> is 
/(x, p, f) - F{xo, po) where xo, po are the phase-space coordi- 
nates at f = of the characteristics that at time t reach x, p 
and F denotes the assigned initial condition for the (positive- 
definite) distribution function at f = 0, i.e. /(xq, po, f = 0) = 
F(xo, Po). Note that dx dp - dxo dpo, i.e., the Jacobian of the 
coordinate transformation in phase space from x, p to Xq, po is 
equal to unity (corresponding to the invariance properties of 
one of the so-called Poincare integrals). Clearly in the xq, po 
coordinates the Vlasov equation reads 



dt 



= 0, 



(A.3) 



where the partial time derivative is taken at constant Xq, po- 

Suppose we are interested in studying the non-conservation 
of a given quantity Q (which might be the single-star energy or 
radial action). Then we may introduce a different coordinate 
system in phase space where instead of X(),po we introduce 
the pair go, 2, where Q = Q{t) = Q{x{t), p(f), is the ap- 
propriate function of the phase space coordinates x, p and of 
time and Qq - Q{t - 0), supplemented by four initial coordi- 
nates, which we symbolically denote by aj, with j - 1,..,4. 
If Q is not conserved, we expect that for t the Jacobian 
J' = STiQo, Q, ctf, Xo, Po) be different from zero; but at f = 
we have Q - Qo, so that the coordinate transformation be- 
comes singular and J' - Q). Clearly, if Q were a conserved 
quantity the above coordinate transformation would be singu- 
lar at all times; a similar situation may occur in the presence of 
an adiabatic invariant, when a quantity Q, such as the single- 
particle energy, is uniquely determined once Qq, aj, and f are 
given. In addition, the Jacobian J' may vanish at special values 
of f, for special values of Qo, Q and aj, when a "caustic" type 
behaviour occurs. 

In terms of the coordinates Qo, Q, aj the Vlasov equation 
for F = F(Qo, Q, aj, t) = F{xo, po) reads 



OF ^ DQ(t) d^^Q 
dt ^ Dt dQ 



(A.4) 



where DQ(t)/Dt is to be expressed as a function of Qo, Q, aj, 
and f. Here the partial time-derivative is taken at constant 
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Qo, Q, Qj. In terms of F, the mass density /(x, p, t) dx dp be- 
comes 

f(x,p,t)dxdp = -dQodQdaj. (A.5) 

JiQo, Q,aj) 

The general solution of Eq. iAAl can be written as 

F(Qo, Q, aj, t) = r{QQ, aj, K,^(Q^u Q, O), (A.6) 



where the characteristic Ka (Qo, Q, t) is a solution of Eq. \AA\ 
at fixed Oj. 

An explicit solution of the above equations is given at the 
end of this Appendix for the special case of a one-dimensional 
motion in a time dependent but spatially uniform force. Here 
we consider two paradigmatic examples. Let us first consider 
the case of 



Dm ^ {Q - Qo) 

Dt t 



(A.7) 



which is expected to hold as typical behavior in the initial 
stages of the evolution, i.e., for f — > 0. At fixed Qq, the charac- 
teristic is then given by 



KaSQo,Q,t) 



{Q - Go) 



(A.8) 



so that the general solution for F should be written as a function 
T of the variables Qo, aj, and Ka/. 



HQo, Q, aj, t) = t{Qo, aj, (Q - Qo)/t). 



(A.9) 



Equation ( IA.9> provides the generic form of the distribution 
function in evolution variables for early times such that (Q - 
Qo)IQo 1 and shows that if the initial distribution function 
is only a function of Qo and of aj, it remains independent of Q 
at all times (over the allowed domain of Q at time t). 
In the second example we take 



DQ(t) ^Dln{aj{Qo,aj,t)) 



Dt 



we have 



Dt 



Ka,(Qo,Q,t)^Q/c 



(A. 10) 



(A.ll) 



where oj - coiQo, aj, t) is an effective "frequency" and the al- 
lowed range of Q at fixed Qo, aj is restricted to the single value 
Q = ojQo/coo, with a»o = to{Qo,aj,t - 0). Therefore, in this 
case we find Q - Qo - QoipJlc^o - 1), which compared with 
Eq. jA.9> . illustrates the fact that in the presence of an adiabatic 
invariant the transformation from x, p, t to Qo, Q, aj, t is singu- 
lar at all times. In turn, finding a high concentration of points 
at all times along a curve in a scatter plot would suggest that 
the transformation be close to singular, which provides us with 
a global test of adiabatic invariance. 

The scatter plots mentioned in Subsection l2.4l and used in 
Sect.l^thus correspond to the mass density D{Qo, Q, t) in the 



{Qo, Q) plane obtained by integrating at time t the mass density 
in phase space over the coordinates Uj 

r{Qo,aj,K„^{QQ,Q,t)) 



D(Qi 



),Q,t)- J doj 



J(Qo,Q,aj) 



For f ^ we have 



D(Qo,Q,t)^ 



!D(Qo,(Q-Qo)/t) 



(A. 12) 



(A. 13) 



The framework that we have developed suggests a number 
of tests and of scatter plots based on different variables, espe- 
cially in the limit f — > 0. We postpone these investigations to 
a separate paper, because they touch on issues that are not of 
immediate astrophysical interest, beyond the main goals of this 
article. 

To conclude, we exemplify the theoretical framework just 
outlined by considering one simple time-dependent system for 
which the characteristics can be easily computed. 

We consider the one-dimensional motion of a star in a ho- 
mogeneous time dependent gravitational field described by the 
Hamiltonian 

Hix, p, t) ^ p^l2- 2xt. (A. 14) 

The characteristics are easily computed and are given by 



P = Po + r, 

X — xo + Pot + /3, 

so that 



H(xo, po, t) - Pq/2- 2xot - pot^ - 



(A. 15) 



(A. 16) 



Note that Ho = Hoixo, po, f = 0) = pl/2. 

We now focus on the non-conserved quantity H (which was 
called Q above) and consider the transformation (xo, po) — » 
(Ho, H), which is easily recovered from: 



PO = {2Ho)'l\ 

xo = (Ho - H)/(2t) - i2Ho)"^t/2 



t^/12. 



(A. 17) 



Here the square roots are properly defined so as to keep track 
of the sign of po- From Eq. (IA.16> we obtain 



DH 2t^ dH 

— ^-2xo-2pot--^-2x^—, 

which can be rewritten as 

DH _H-Ho 1/2 P 



(A.18) 



(A. 19) 



and reduces to DH/Dt ~ {H - Ho) It for f 0. The character- 
istic function K(Ho, H, t) defined by Eq. iAAl is given by 

K{Ho,H,t) = (Ho- H)/(2t)-(2Hoy'^t/2-t^/l2 

= xoiHo,H,t) (A.20) 
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as can be easily verified by substitution and is a priori obvi- 
ous from the fact that, at constant Hq, the characteristics are 
labelled by the value of xq. 

The Jacobian J'iHo, H) is given by 



.20/2 



(A.21) 



Thus from the distribution function F(xo,po), in the H,Ho,t 
variables we have 



F{xo,po) dxodpo = 



r(K(Ho,HM2Ho)"^) 



dHodH, 



(A.22) 



with K(Ho, H, t) defined by Eq. jA.20> . Note that, if at f = the 
distribution function is homogeneous in space it is constant in 
H at all times. This is best realized by considering a distribution 
that is Gaussian both in po and in xo, for which we obtain in 
H, Ho, t 
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which becomes independent of H in the limit L — > 00. For 
finite L, in the limit f — > 00 the distribution function is centered 
around H = and has a width that increases much slower, 

as f i.e., the distribution tends to become "monochromatic". 

Two different examples, the case of a parametric linear os- 
cillator and the case of a free star in a small- amplitude wave- 
like potential, have also been solved explicitly in this frame- 
work, demonstrating how adiabatic invariance, the formation 
of caustics, and the effects of resonances may be recognized by 
means of studies of the Vlasov equation in evolution variables. 
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Fig. 1. Case A (dynamical friction): collisional shell of Nf frag- 
ments inan« = 3 polytropic galaxy (runAl). The top panel de- 
scribes the time evolution of the radii of the spheres containing 
15, 35, 55, 75, 95 % of the total shell mass. The central panel 
compares the final galaxy density profile to that of the corre- 
sponding control case, which is very close to the initial galaxy 
density profile. The bottom panel illustrates the time evolution 
of the galaxy global anisotropy parameter {IKrlKj). 



Fig. 2. Cases B (slow and adiabatic) and C (fast, non-adiabatic): 
evolution of an n = 3 polytropic galaxy in the presence of an 
external potential which mimics the fall of the shell observed 
in run Al. The three panels describe the same quantities as in 
Fig. [2 Note the qualitative difference with case A in the evo- 
lution of the density profile. In the middle frame the B and C 
curves are basically identical. In the bottom frame the fast case 
CI is illustrated only in the insert. 
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Fig. 5. Comparison of final total density profiles in differ- 
ent evolution processes (solid lines, dynamical friction; dot- 
ted lines, control cases; dash-dotted lines, adiabatic evolution) 
for the polytropic galaxy model (upper frame) and for the f-'^^ 
model (lower frame). The control cases, which are basically 
time-independent, thus illustrate the adopted initial conditions. 



